Magnetometry via a double-pass continuous quantum measurement of atomic spin 
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We argue that it is possible in principle to reduce the uncertainty of an atomic magnetometer 
by double-passing a far-detuned laser field through the atomic sample as it undergoes Larmor pre- 
cession. Numerical simulations of the quantum Fisher information suggest that, despite the lack of 
explicit multi-body coupling terms in the system's magnetic Hamiltonian, the parameter estimation 
uncertainty in such a physical setup scales better than the conventional Heisenberg uncertainty 
limit over a specified but arbitrary range of particle number N . Using the methods of quantum 
stochastic calculus and filtering theory, we demonstrate numerically an explicit parameter estimator 
(called a quantum particle filter) whose observed scaling follows that of our calculated quantum 
Fisher information. Moreover, the quantum particle filter quantitatively surpasses the uncertainty 
limit calculated from the quantum Cramer-Rao inequality based on a magnetic coupling Hamilto- 
nian with only single-body operators. We also show that a quantum Kalman filter is insufficient to 
obtain super-Heisenberg scaling, and present evidence that such scaling necessitates going beyond 
the manifold of Gaussian atomic states. 
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I. INTRODUCTION 

It is well-appreciated in physics that the properties of 
a field must often be determined indirectly, such as by 
observing the effect of the field on a test particle. Take 
magnetometry for example: the strength of a magnetic 
field might be inferred by observing Larmor precession 
in a spin-polarized atomic sample P and estimating the 
field strength B from the precession rate. Inherent in this 
process is the fact that the atomic spin must be measured 
to determine the extent of the magnetically-induced dy- 
namics. For very precise measurements, uncertainty 5B 
in the estimated value B of the field is dominated by 
quantum fluctuations in the observations performed on 
the atomic sample. The results presented here fall under 
the umbrella of quantum parameter estimation theory 
m [3], where the objective is to work within the rules 
of quantum mechanics to minimize, as much as possible, 
the propagation of this quantum uncertainty into the de- 
termination of metrological quantities, like B. 

Given, for instance, a ?/-axis magnetic field B = 5y, 
an atomic sample couples to B via the magnetic dipole 
Hamiltonian 



H 



-hjBFy 



(1) 



where 7 is the atomic gyromagnetic ratio and Fi — 

X^jii /P^ = ^7 y? 2) the collective spin operators 
obtained from a symmetric sum over N a identical spin- 
/ atoms. If the atoms are initially polarized along the 
X-axis, the Larmor dynamics and thus B can be in- 
ferred by observing the z-component of the atomic spin 
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Fz [iiHiiaE]. 

It is possible to place an information-theoretic lower 
bound on the units-corrected mean-square deviation of 
any magnetic field estimator B from the actual value of 
the field B, 



SB 



B 



\\d{B)/dB\ 




(2) 



via the quantum Cramer-Rao inequality (as described in 
Section |lll| H H [3 [8]. The behavior of the estimator 
uncertainty with the number of atoms Na depends on 
the characteristics (e.g., separable, entangled, etc.) of the 
quantum states used to compute the expectation value in 
Eq. ([2| as well as the nature of the induced dynamics [9] . 
If one does not permit quantum entanglement between 
the different atoms in the probe, it can be shown that 
the optimal parameter resolution obtained from Eq. ([T]), 
evolved for a time t, is given by the so-called shotnoise 
uncertainty [lj|5] 



(3) 



whose characteristic 1/\/Na scaling is a byproduct of 
the projection noise (AFz) = \/?y2 for a spin coherent 
state [To] (here F = fNA for a sample of Na atoms each 
with total spin quantum number /). It was believed for 
some time that the fundamental limit to parameter es- 
timation, even when exploiting arbitrary entanglement 
between atoms in the probe, offers only a quadratic im- 
provement 



a 



(4) 



up to an implementation-dependent constant a. Eq. Q 
has traditionally been called the Heisenberg uncertainty 
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scaling, and it can be achieved in principle for various 
spin resonance metrology problems ^10^, including mag- 
netometry [5]. For an ensemble of A^^ spin- 1/2 particles 
prepared into the initial cat-state (|tT " • • T) + |ii • • • i 
))/v/2, the uncertainty scaling is given by I/^INa and is 
sometimes called the Heisenberg Limit. 

Recently, however, it was shown that 1/Na scaling can 
be surpassed! 9 by extending the linear coupling that un- 
derlies Eq. ([ij to allow for multi-body collective interac- 
tions [9l [11] . Were one to engineer a probe Hamiltonian 
where B multiplies /c-body probe operators, such as Fy , 
then the quantum Cramer-Rao bound [3 indicates that 
the optimal estimation uncertainty would scale more fa- 
vorably as AB]^ ~ i^.' Unfortunately, metrological 
coupling Hamiltonians are rarely up to us — they come 
from nature, like the Zeeman interaction — suggesting 
that one is stuck with a given uncertainty scaling with- 
out changing the fundamental structure of Eq. ([T]) . Fur- 
thermore, it was shown in Ref. ^ that the addition of 
an auxiliary parameter- independent Hamiltonian Hi{t) 
such that 



H = -h-fBFy^Hi{t) 



(5) 



does not change the scaling of the parameter uncertainty 
for any choice of Hi{t). 

At the same time, however, it should be well- 
appreciated that the dynamics one encounters in any ac- 
tual physical setting are effective dynamics. Indeed, even 
the hyperfine Zeeman Hamiltonian Eq. ([T]) is an effective 
description at some level. This begs the question as to 
whether one can utilize an auxiliary system to induce 
effective dynamics that improve the uncertainty scaling 
in quantum parameter estimation by going outside the 
structure of Eq. ([5|. The purpose of this paper is to 
provide some direct evidence that doing so is possible. 

We study effective nonlinear couplings generated by 
double-passing an optical field through an atomic sam- 
ple (q.v. Figure [T]) [12^ JHl [M]. Continuous measure- 
ment of the scattered field then allows for the estimation 
of Fz and by extension, the magnetic field. Building on 
the quantum stochastic calculus approach in [l^, we de- 
rive the quantum filtering equations for estimating the 
state of the atomic sample. Although the effective dy- 
namics are no longer described by a Hamiltonian, we 
perform numerical calculations of the quantum Fisher 
information to obtain a theoretical lower bound on the 
uncertainty scaling of an optimal magnetic field estimator 
[3\ Our simulations suggest that for certain parameter 
regimes, the double-pass system's sensitivity to magnetic 
fields scales better than that of a comparable single-pass 
system and what would be computed by applying the 
methods of Ref. [9 to Eq. ([5|. We also conduct direct 
simulations of magnetic field estimation for our system 
using a general nonlinear parameter estimator, called a 
quantum particle filter [15], as further evidence for the 
improved uncertainty scaling provided by our proposed 
magnetometer. 



Unfortunately our results are somewhat muted by the 
fact that we have, despite our best efforts, not yet 
found a parameter estimator whose uncertainty scaling 
can be shown analytically to outperform the conven- 
tional Heisenberg limit. In particular, we show that im- 
proved scaling is not achieved by a quantum Kalman fil- 
ter [ini [171 118] , as such a filter is only suitable for estimat- 
ing magnetic fields in the linear small-angle regime and 
where the state is Gaussian and the dynamics are well 
approximated by a low order Holstein-Primakoff expan- 
sion [51 [19]. Although Kalman filters have had success 
in describing the single-pass system [5 , simulations sug- 
gest the Gaussian and small-angle approximations break 
down precisely when exact simulations of the double-pass 
system show improved sensitivity. For pedagogical pur- 
poses, we detail the derivation of such linear- Gaussian 
filters using the method of projection filtering [21]. 
Doing so allows us to observe directly the limitations 
that arise when imposing the small-angle and Gaussian 
assumptions, and it also provides a framework for the 
future development of more sophisticated filters. As a 
consequence, our analysis here is restricted to values of 
the total atomic angular momentum F that is small com- 
pared to experimentally relevant values. It remains an 
important problem to develop nonlinear parameter esti- 
mation methods suitable for implementation in an actual 
laboratory setting. 

The remainder of the paper is organized as follows. 
In Section [Hj we detail the magnetometer setup and its 
corresponding quantum stochastic description. We also 
present the quantum filtering equation which describes 
continuous measurement of the double-passed atomic 
sample. Using this model, we present a numerical anal- 
ysis of the quantum Cramer- Rao bound in Section HI) 



and observe an improved uncertainty scaling relative to 
the shotnoise and Heisenberg scalings. In Section [IV| we 
use the double-pass quantum filter to develop quantum 
particle filters and quantum Kalman filters suitable for 
estimating an unknown magnetic field. In Section [V| we 
present simulations of the particle filter magnetic field es- 
timator and discuss the evidence for improved magnetic 
field sensitivity in the double-pass system. We conclude 
in Section [Vll 



II. 



CONTINUOUS MEASUREMENT OF THE 
DOUBLE-PASS SYSTEM 



Consider the schematic in Fig. [T] The objective of this 
apparatus is to estimate the strength of a magnetic field 
oriented along the laboratory ?/-axis by observing the ef- 
fect of that field on the spin state of the atomic sam- 
ple. Like most atomic magnetometer configurations, our 
procedure relies upon Larmor precession and uses a far- 
detuned laser probe to observe the spin angular momen- 
tum of the atomic sample. Unlike conventional atomic 
magnetometer configurations, however, the probe laser 
is routed in such a way that it passes through the atomic 
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the familiar Stokes operators: 



FIG. 1: (color online) Schematic of a broadband atomic mag- 
netometer based on continuous observation of a polarized op- 
tical probe field double-passed through the atomic sample. 



sample twice prior to detection [T2l [131 E] • 

Qualitatively, the magnetometer operates as follows. 
The incoming probe field propagates initially along the 
atomic z-axis and is linearly polarized. As a result of 
the atomic polarizability of the atoms, the probe laser 
polarization acquires a Faraday rotation proportional to 
the z-component of the collective atomic spin. Two fold- 
ing mirrors are then used to direct the forward scattered 
probe field to pass through the atomic sample a second 
time, now propagating along the atomic ?/-axis. Prior to 
its second interaction with the atoms, polarization op- 
tics convert the initial Faraday rotation into ellipticity. 
Thus on the second pass, the atoms perceive the opti- 
cal helicity as a fictitious ?/-axis magnetic field acting in 
addition to the real field 5, providing a positive feed- 
back effect modulated by the strength of B. That is, the 
total Larmor precession of the spins is enhanced by an 
amount which depends on the true strength of B. The 
twice forward-scattered optical field is then measured in 
such a way that is sensitive only to the Faraday rotation 
induced by the first pass atom-field interaction. 



A. Quantum Stochastic Model 



When the collective spin angular momentum of a mul- 
tilevel atomic system interacts dispersively with a trav- 
eling wave laser field with wavevector k, the atomic spin 
couples to the two polarization modes of the electromag- 
netic field transverse to k. These polarization modes can 
can be viewed as a Schwinger-Bose field that when quan- 
tized in terms of a plane- wave mode decomposition yields 



1 /.f ^ ^ 

2 y^-\-,u^-\-,uj ~^ ^-,u^-,uj 



(6) 



2 



(8) 



(9) 



Here, we have expressed the Stokes operators in terms of 
the Schrodinger-picture field annihilation operators, ^ 
and «y^^, for the plane- wave modes with frequency co and 
linear polarization along the x- and y-axes, respectively, 
as well as their corresponding transformations into the 
spherical polarization basis. 

In developing a physical model for the atom-field in- 
teraction in Fig. [l] it is convenient to transform from 
a plane-wave mode decomposition of the electromagetic 
field to operators that are labeled by time. Towards this 
end, we define the time-domain Schwinger boson annihi- 
lation operator as the operator distribution 



st = - g{u^) al^^a^^y'^duj, 

^ J-oo 



(10) 



where g{u;) is a form factor. This definition permits us 
to express the Stokes operators as 

^z,t = i (4 - and Sy^t = - (^t + 4) ^ (H) 

which should be reminiscent of quadrature operators and 
also places the field operators in a form that is directly in 
line with the standard nomenclature adopted in the field 
of quantum stochastic calculus. 

With a suitable orientation of the polarization optics 
(A/2 and A/4) in Fig. [l] the interaction Hamiltonians for 
each pass of the probe light through the sample are then 

= +n/iF,5,,t = +^^M^z(4-^t) (12) 



hnP^ (st + 4) , (13) 



respectively, where the coupling strengths /i, n arise from 
the atomic-polarizibility level structure for the particular 
atoms being used [22> ,23^ . Note that in developing these 
Hamiltonians, which are of the standard atomic polar- 
izability form, it was assumed the laser frequency uoi is 
far detuned from any atomic resonance so that rank-two 
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spherical tensor interactions can be neglected [22j |23] . In 
practice, the validity of such an assumption can depend 
heavily on the choice of atomic level structure and exper- 
imental parameters such as the intensity and detuning of 
the probe laser field. 

In addition to specifying the Hamiltonians for the two 
atom-field interactions, it is also necessary to stipulate 
the measurement to be performed on the probe laser. 
Since we expect that the amount of Larmor precession 
(possibly augmented by the addition of the double-passed 
probe field) will cary information about the magnetic 
field strength 5, we must choose the measured field op- 
erator Zi appropriately. Since the magnetic field drives 
rotations about the atomic y-axis, it is the 2;-component 
of the atomic spin that indicate such a rotation. From 
the form of the first-pass interaction Hamiltonian il^^ ^ 
we see that the z-component of the atomic spin couples to 
dynamics generated by the field operator s^^i = i[s\ — s^). 
The affect of such a coupling is then observed by mea- 
suring the orthogonal quadrature, indicating that the ap- 
propriate polarization measurement should be Sy^f 



1. The Stochastic Propagator 

Analyzing the two individual interactions Hi and H2 is 
a well studied problem in continuous measurement theory 
as a weak-coupling limit [24^ 25, 26, 27 . Heuristically, 
we are interested in the long time cumulative effects of 
the fast field dynamics relative to the slow atomic dy- 
namics. Much like the classical central limit theorem, 
the above references show that the cumulative effect of 
the field converges to a white noise process, so that the 
overall interaction-picture evolution of the joint atom- 
field system is well-described by the following quantum 
stochastic differential equation (QSDE) propagators 

dU^^^ = ^^^MFz{dSl - dSt) - ^mP^dt - ^Hdt"^ U^^^ 

(14) 

dU^^^ = i^i^/kFy{dSt + dSl) - ^kP^dt - ^Hdt"^ /7f ^ 

(15) 

As discussed in ^24j, the rates m = 27r|/i^(a;/)p and k = 
27r\ Kg{u;i)\'^ are evaluated in terms of the original coupling 
strengths and form factor at the drive laser frequency 
uJi^ H is an arbitrary atomic Hamiltonian and dS} and 
dSt are delta-correlated noise operators derived from the 
quantum Brownian motion 



f 

Jo 



,,du. 



(16) 



The noise terms satisfy the quantum Ito rules: dStdS] = 
dt and dSjdSt = dS^ = {dS})'^ = 0, and can be viewed 
heuristically as a consequence of vacuum fiuctuations in 



the probe field. Note that these are field operators and 
are defined so as to commute with all operators Xs for 
s<t. 

To obtain a single weak-coupling limit for the double- 
pass interaction, we combine the separate equations of 
motion for the two propagators into a single weak- 
couping limit as follows. First, write the two single-pass 
evolutions in terms of the generators of the dynamics 



da 



(1) 



dtU^ 



(1) 



and, dUi 



(2) 



(1) 



(17) 



and then expand the differential dUt of the combined 
propagator 

dUt+6t = {i^bt){i^dt)Ut (18) 
= Ut^ [d^b^bd^ Ut (19) 

such that the combined propagator dUt = Ut-\-st — Ut then 
satisfies 



dUt 



(d^b^ bd^ Ut 



(20) 



After evaluating the combined evolution for the propa- 
gators in Eqs. (14) and (Il5| in light of the quantum Ito 



rules, we find that the single weak-coupling limit propa- 
gator satisfies 



ivkmFyFzdt ■ 



1 



mP^dt 



-kP^dt 
2 ^ 



2i 



-Hdt 



mP^{dSl - dSt) + iVkPy{dSl + dSt)] Ut^ (21) 



dUt 



Observe that as a result of the manner in which the 
combined weak-coupling limit was taken, the Hamilto- 
nian term has the property that rates which appear in 
it differ by a factor of two from those that would be ex- 
pected from a single weak-coupling limit. This factor of 
two is essentially the rescaling of time units that arises 
from aggregating two sequential weak-coupling limits as 
a single differential process. To retain consistency with 
the original definition of the frequencies that appear in 
the parameter-coupling Hamiltonian, it is essential to 
rescale time units such that frequencies in the parameter- 
coupling Hamiltonian are as expected. Doing so is ac- 
complished by reversing the effective 2dt dt transfor- 
mation that occurred in the derivation, and thus dividing 
all rates by two to give 



dUt 



iVKMpyP^dt - -Mp^dt - -Kp^dt - iHdt 
2 2 ^ a 



Mp^{dS\ - dSt) + iVKPyidSl + dSt) Ut (22) 

where M = m/2 and K = k/2. We note that this final re- 
sult agrees with the propagator obtained by Sarma et. al 
[13 who also derived the quantum stochastic propagator 
of this system in order to characterize the generation of 
polarization and spin squeezing as suggested by Sherson 
and M0lmer [12^ and Muschik et. al [14 . 
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B. Double-Pass Quantum Filter 

Continuous polarimetry measurement of the transmit- 
ted probe field is well described by a quantum stochas- 
tic calculus model [SQ which in turn allows for statis- 
tical inference or quantum filtering of the state of the 
atomic sample [27l |29] . We now derive the quantum fil- 
tering equation corresponding to the propagator in Eq. 
( 22 ) and continuous observation of the forward-scattered 



probe polarization. Our derivation follows the general 
approach developed in Section 4 of Ref. [27 . In this con- 
densed presentation, we gloss over some of the technical 
details that one would consider in a completely rigorous 
derivation of a quantum filter. In particular, we do not 
demonstrate that our result is optimal, nor even delve 
into the details of what optimal means in this context. 
Suffice it to say, the results hold up under more rigorous 
mathematical scrutiny, and the interested reader is en- 
couraged to consult Refs. [271 [12] for a more thorough 
discussion of these details. Readers who are familiar 
with mathematical filtering theory will miss little by only 
skimming the derivation as it is presented here. 

The task of the quantum filter is to construct a best- 
estimate of atomic observables conditioned on continuous 
observation of the transmitted probe field. Our first step 
is one of convenience — we rewrite the double-pass prop- 
agator of Eq. ( 22 ) in the general form [29 



dUt 



LdSj - PdSf - -PLdt - iHrdt 



Uu (23) 



identifying L = VMF^ + iVKFy and He = H - 
VKM{F^Fy + FyF^)/2. Then, the Heisenberg-picture 
evolution of any atomic observable X is given by the so- 
called quantum flow 



3t{X) = UKX)Uf 



(24) 



Applying the Ito product rule (discussed in Appendix 
[A[ ) twice and noting the Schrodinger picture operator sat- 
isfles dX = 0, we find 



djt{X) = UlXdUt^dUlXUt^dUlXdUt 



(25) 



= j,([X,L])d^J + j,([Lt,X])d^, + j,(£(X))dt, 
where we have defined the familiar Lindblad generator 

C{X) = i[H,X]^ PXL- ^{PlX ^xPl). (26) 

Similarly, the observation process corresponding to po- 
larimetric detection of the ^/-component of the vector 
Stokes operator is given by 



(27) 



Two applications of the Ito product rule, which generate 
seven initial Ito products, gives its time evolution as 



dZt = jt{L + P)dt + dSl + dSt. 



(28) 



In the language of classical control theory, Eqs. (25) 
and ([28]) are a system/observation pair. We see that 
the observations process carries information about the 
state of the system, jt{L + 1/^), albeit corrupted by the 
quantum white noise term dSl-\-dSt. Moreover, as seen in 
Eq. (25), the system itself is also driven by the quantum 
white noises. The quantum filter is tasked with picking 
out the relevant information from the observations and 
combining it with the model of the dynamics in order to 
estimate system observables. 

Mathematically, the filter is given by the conditional 
expectation 

(29) 



7rt[X]=nMX)\Z^o,t)], 



where ^(o,t) is the entire measurement record up to time 
t. Intuitively, we expect that the filter will be some func- 
tion of the measurement process on which we are condi- 
tioning. This statement can be made more mathemati- 
cally rigorous by showing that the conditional expecta- 
tion in Eq. (29) is adapted to the measurement process. 



Furthermore, it is possible to derive a recursive form for 
the filter, in which only the most recent measurement 
increment is needed to update our estimate, and the dy- 
namical equation for the filter can be written as 



d7Tt[X] = kt{X)dt + mt{X)dZt, 



(30) 



where kt{X) and mt{X) are yet to be determined func- 
tions and the form is the most general we could consider, 
given that we only have access to the measurements. In 
order to find the forms of these functions, we can utilize 
the property of conditional expectation 



E[fE[X|Z(o,t)]]=E[fX], 



(31) 



which is true for any Y that is a function of the measure- 
ment process. 



One way of enforcing Eq. (31) is by ensuring that the 
following holds: 



E TTt 



-t[X]efo9is)dz. 



= ]E\jt{X)e^o9is)dz^ ^ ^32) 



since by taking linear combinations and appropriate 
derivatives of either side, any (analytic) function Y of 
the measurement current can be generated via an appro- 
priate choice of g{t). To simplify the algebraic manipu- 
lations, we multiply both sides by the deterministic inte- 
grating factor exp(— ^ g{s)'^ds) and identify our "gen- 
erating function" as 



exp 



g{s)dZs - \ f g{sfds ] . (33) 
^ Jo 



To calculate the stochastic differential equation for ef, 
we find it useful to set dRt = g{t)dZt — ^g{t)'^dt, so that 
ef = ex.p{Rt). Using the chain rule and expanding to 
second order in differentials, we have 

1 



del 



el dRt - 
g{t)eldZt. 



-elidRtY 



(34) 
(35) 
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where we have used the fact that (dRt)'^ = {g{t)dZt)'^ = 
g{t)^dStdSl = g{t)^dt 

Since the individual white noise terms, dSt^dS}, are 
zero in expectation, we can formahy take the time deriva- 
tive of either side of Eg. (32) as long as we use the Ito 



product rule of Eq. ( A3 ) prior to taking the expectation. 
Doing so, we find 



E 



+g{t)el (mt{X)+jt{L + L^)7rt[X] 



dt 



nelMX)] = E eljt{C[X]) 



(37) 



+g{t)e!jt{L^X + XL) 



Note that in reaching the above form of Eq. (37), we 
made use of the identities 



jt{A)jt{B) = MAB) 
jt{A)-jt{B)=jt{A-B), 



(38) 
(39) 



which are readily verified using the definition of jt{X). 

To obtain the quantum filtering equation, we equate 
like terms in Eqs. (36) and (l37|) to enforce Eq. (I32|), and 



to ensure that rrit and kt remain functions of the mea- 
surement record, we replace ]E[j^(-)] with ]E[7r^(-)], as this 
is just Eq. (32) with g{t) = 0. After a little rearranging. 



we arrive at the quantum filtering equation 

dntiX] = nt[C[X]]dt (40) 
+ (ntiPx + XL] - TTtiP + L]7rt[X]) 

X (dZt - TTtlL ^ P]dt^ 

Oftentimes, and particularly in the traditional quan- 
tum optics setting, one works with the so-called adjoint 
form for the filter, which gives a dynamical equation for a 
density matrix pt that satisfies 7rt[X] = Tr [Xp^] . This is 
easily determined from Eq. (|40|. Plugging in the specific 
L and H for our double-pass system undergoing Larmor 
precession, we find the double-pass quantum filter 



dpt = i-fB[Fy,pt]dt^iVKM[Fy,{F,,pt}]dt 

^MV[F,]ptdt + KV[Fy]ptdt (41) 
+ ^VMM[F,]pt^iVK[Fy,pt]^ dWt 

where the innovations process 



dWt = dZt - 2VMTT[F^pt]dt 



(42) 



is a Wiener process, i.e. E[(ilVt] = 0,dW^ = dt. The 
various superoperators are defined as 

Vmpt = F.ptPl - ^FlFkPt - \ptFlFu (43) 

M[F,]pt = F,pt + ptF, -2Ti[F,pt\pt (44) 

{F,,pt} = F^Pt^PtF^ (45) 



In looking at the structure of the filter, we see a deter- 
ministic term which evolves the current state estimate via 
the standard open system evolution one would expect for 
the system. The stochastic term performs the condition- 
ing and is weighted by the innovations process, which 
captures how much new information the latest measure- 
ments provide. 

One other form which is useful when the quantum 
state remains pure is the stochastic Schrodinger equa- 
tion (SSE). As developed in Appendix [bJ the SSE for 
the double-pass quantum filter is 



U 



ijBFy 



M 



(F. - {F^)tf 



\''~~' 2 
-^iVKMFy{F,^{F, 

+ (Vm(F, - {F,)t) + iVKFy^ \^)tdWt. 



(46) 



] Wtdt 



III. 



THE QUANTUM CRAMER-RAO 
INEQUALITY 



In order to characterize the performance of the mag- 
netometer, we consider quantum information theoretic 
bounds on the units-corrected mean-square deviation of 
the magnetic field estimate B of the true magnetic field B 
[31 H], given in Eq. (2|. The quantum Cramer- Rao bound 
[11 m d E] states that the deviation of any estimator is 
constrained by 



SB > 



lB{t) = TT[pBmlit)]. (47) 



where the quantum Fisher information Xsit) is the ex- 
pectation of the square of the symmetric logarithmic 
derivative operator, defined implicitly as 

= \{ZB{t)pB{t) + PBmB{t)). (48) 

For pure states, p\ = pB-, so that 



S.B{t) = 2- 



dB 



which indicates 



5B > 



1 (dpBjt) 

2 \ V dB 



2\ ~2 



(49) 



(50) 



In this form, we see that the lower bound is related to 
the sensitivity of the evolved state to the magnetic field 
parameter. That is, any estimator's performance is con- 
strained by how well the dynamics transform differences 
in the value of B into differences in Hilbert space. 

As discussed by Boixo et. al in [9 , for Hamiltonian 
evolution, the quantum Cramer-Rao bound may be ex- 
pressed in terms of the operator semi-norm, which is 
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the difference between the largest and smahest (non- 
degenerate) eigenvalues of the probe Hamiltonian. For 
the magnetic dipole Hamiltonian in Eq. ([T]), this bound 
is simply the Heisenberg limit in Eq. Q. More gener- 
ally, the authors show that a probe Hamiltonian which 
involves /c-body operators gives rise to an uncertainty 
scaling of 1/tF^. They further argue that no ancillary 
quantum systems or auxiliary Hamiltonians contribute 
to this bound; it is determined solely by the Hamiltonian 
that directly involves the parameter of interest. 

Such analysis suggests the double-pass quantum sys- 
tem, whose only direct magnetic field coupling is in the 
magnetic dipole Hamiltonian, should show no more sen- 
sitivity than a single pass system. There are several rea- 
sons why we believe there is more to the story. Firstly, 
the unitary evolution of the joint atom-field system in 
Eq. (23) involves an auxiliary system of infinite dimen- 



sion. As such, it is not clear that the arguments leading 
to the operator semi-norm are valid, in particular due 
to the fact that the white noise terms dSt , dS} are sin- 
gular. Additionally, the double-pass limit is a Markov 
one, in which the interaction the light field mediates be- 
tween atoms is essentially instantaneous relative to other 
time-scales in the problem. The effective interaction is 
therefore fundamentally different than one in which mea- 
surements of a finite dimensional ancilla system are used 
to modulate the evolution of the probe atoms. Thus the 
conditioned system, given in terms of the quantum fil- 
ter of Eq. (40) or (46), does not correspond to unitary 



dynamics. Indeed, looking at Eq. (46), we see that the 
local generator of dynamics is path- dependent^ given in 
terms of the expectation of Fr^. Therefore, as the mag- 
netic field directly impacts the state through the mag- 
netic dipole term, it also non-trivially modulates future 
dynamics through a state-dependent generator. 



A. Numerical Analysis of the Quantum Fisher 
Information 



Unfortunately, it is not clear how to fold the quan- 
tum stochastic or quantum filtered dynamics analytically 
into the semi- norm bound considered in 9 . Nonethe- 



less, the quantum Cramer- Rao bound in Eq. (50) is ex- 



cellent fodder for computer simulation. By numerically 
integrating the stochastic Schrodinger form of the quan- 
tum filter in Eq. (46), a finite difference approximation of 



dpB{t)/dB may be evaluated for different collective spin 
sizes F. That is, for a given choice of F, a finite difference 
approximation of the quantum Fisher information near 
B = can be constructed by co-evolving three trajecto- 
ries, po{ts)^ p5B{ts)^ and p -5 Bits) (seeded by the same 
noise realization), and calculating 




p5B{ts) - p-5B{ts) 

2SB 



The expectation value is then evaluated numerically by 
averaging over a statistical ensemble of such measure- 
ment trajectories. 



1. Simulation Results 

Figure |2] shows the results of such simulations of the 
quantum Fisher information as a function of the spin size 
F and different relative values of the coupling strengths 
M and K. For convenience, all rates are taken to be unit- 
less, defined in units of an arbitrary frequency v (which is 
then set to = 1 in the actual simulation). For all cases 
described here, the quantum filter in Eq. (46) was inte- 



Po(^s)]. 

(51) 



grated from an initially x-polarized spin coherent state 
using a second-order predictor-corrector stochastic inte- 
grator [30 with a step size dt = 10~^z^~^ until a time 
ts = 0.1h'~^. The density operators in Eq. (51) were con- 
structed as pB{ts) = \i^B{ts)){i^B{ts)\ for the state of the 
system at time ts evolved under the magnetic field with 
value B. The finite-difference parameter was chosen to be 
SB = 5 X 10~'^u and 7 = 1. Due to the stochastic nature 
of the evolution, the finite difference approximation was 
averaged over 100 trajectories, with each member of the 
three density matrices coevolved using the same white 
noise realizations. The plots in Fig. [2] correspond to the 
lower bound determined by Eq. ([50| using the finite dif- 
ference approximation and finite averaging approxima- 
tion just described. The error bars correspond to the 
standard deviation of the lower bound for the ensemble 
of 100 trajectories. A similar finite difference approxi- 
mation around the mean value B = O.lu was conducted 
and showed the same behavior as that described here for 
a mean value of B = 0. 

Figure [2|a) plots the estimator uncertainty lower 
bound as a function of F for a double-pass system, with 
M = lOu and K = 6 x 10~^u^ and an equivalent single- 
pass system with M = lOu and K = 0. For these param- 
eter values, the simulations show a clear improvement 
in the scaling of the uncertainty lower bound, even tak- 
ing into account the increased fiuctuations of the double- 
pass system. Also shown are the analytic scalings for 
shot-noise uncertainty Eq. (|3|, Heisenberg uncertainty 
Eq. (|4| (with a = 1) and the quantum Cramer-Rao 
bound for a hypothetical magnetic field coupling inter- 
action — hjBFy. A least-squares fit to a power-law form 
suggest the single-pass system scales as F~^-^^ ^ which is 
in good agreement with the Heisenburg-limited scaling of 
F~^. A power law fit for the double-pass system gives 
an exponent of —1.34, although its lower bound appears 
to have more complicated scaling. In fact, the behavior 
towards larger F suggests that the double-pass system 
may show a super-polynomial uncertainty scaling. 

We note that there is a constant prefactor difference 
between our estimator simulations and the conventional 
Heisenberg limit. This prefactor appears in the plot as a 
constant offset between the simulated single-pass bound 
and the l/ts2F Heisenberg bound. There are two reasons 



8 




FIG. 2: (color online) Comparison of the lower bound on estimator uncertainty for the single- and double-pass systems 
dete rmined by finite- difference calculations of the quantum Fisher information. Both plots approximate the derivative in Eq. 
(50) using SB = 5 x 10~^iy at a final integration time ts — 0.1z^~^. The lower bound estimate is averaged over 100 simulated 
trajectories with the error bars indicating the standard deviation observed in those 100 samples, (a) Lower bound scaling for 
M = lOz/, = 6 X 10"^i/. A least -squares fit shows that the single-pass uncertainty is consistent with the Heisenberg scaling, 
with a numerical power law fit of , The double-pass clearly improves faster than 1/F over the range of F considered 

here. Also shown for comparison are the analytic shot-noise, Heisenberg and two-body coupling lower bounds, (b) Simulations 
were also performed for a difference choice of double-pass coupling parameters M = K = Iz/, where it can be seen that the 
improved scaling is highly dependent on selecting appropriate coupling strengths. 



why one might expect such an offset: (1) the Heisen- 
berg limit is computed from the Cramer-Rao inequal- 
ity by optimizing over all initial states, including opti- 
mally squeezed states or "cat" superposition states which 
achieve the bound while our procedure begins from a 
spin coherent state which is then gradually entangled by 
the measurement process; (2) a continuous measurement 
scenario based on the theory of open quantum systems 
cannot be expressed in the form of a simple Hamiltonian. 

Furthermore, our simulations reveal that the uncer- 
tainty scaling of the initially x-polarized spin coherent 
state depends sensitively on the particular choice of cou- 
pling strengths M and as seen in Fig. |2|b). This 
figure considers the case of equal coupling strengths, 
M = K = and we see that the double-pass system 
quickly loses its advantage after F ^ 25. This can be un- 
derstood in light of the quantum stochastic model of the 
previous subsection. In considering the general stochastic 
propagator of Eq. ( 23 ) , we identified the coupling opera- 
tor L = y/MFz + iy/KFy, which \i M = K, is essentially 
the angular momentum lowering operator along x — F-^x- 
If M, K 7^, a continuous measurement of this oper- 
ator very quickly moves the +x-polarized initial state 
onto the — x-polarized state, which is an attractive fixed 
point of F-^x' Once this state is reached, the dynamics 
become relatively insensitive to the magnetic field value 
and result in a poor uncertainty lower bound. On the 
other hand, if M, K are much smaller than 7^, the pos- 
itive feedback from the i^/KFy term is washed out by 



Larmor precession. Given that we are interested in de- 
tection limits, i.e. 5 ~ 0, we do not focus on the regime 
where Larmor precession dominates. 

Regardless, Figure [2|b) suggests that for a given M 
and 7^ 0, there will be some value of F after which 
the lowering operator dynamics dominate, rendering the 
double-pass system relatively useless for estimating B. 
But by manipulating M and K relative to F, the extra 
dynamics allow for an improved sensitivity to the mag- 
netic field. Indeed, Fig. |2|a) uses a very small second- 
pass strength relative to the first-pass strength M, so 
that just the right amount of positive feedback enters via 
the i\/KFy term in the measured operator. The positive 
feedback can be viewed as increasing the magnetic fields 
effect on the rate of precession. 

In short, Figure |2|a) suggests that there are some pa- 
rameter values, appropriate for some range of F, which 
show an estimator uncertainty lower bound scaling below 
the Heisenberg limit. It is important to note that com- 
putational constraints limited simulations to F ^ 150, 
which is well below the collective spin size one might ex- 
pect for 10^ — 10^ atoms. However, we believe the scal- 
ings are nonetheless suggestive of quite general behavior 
for this system. Although in practice, it seems that one 
would need to fine tune the coupling strengths M and 
K in order to be in a regime with such scaling. It may 
be that such coupling strengths are inaccessible in an ex- 
perimental setting. While this is an important considera- 
tion, there is a more pressing theoretical question — does 
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a practical estimator exist which saturates the quantum 
Cramer-Rao bound? We summarize our search for such 
an estimator in the fohowing section. 



IV. MAGNETIC FIELD ESTIMATORS 

While studying the properties of lower bounds on esti- 
mator performance is important for developing an under- 
standing of the capabilities of a given parameter coupling 
scheme, any actual procedure for implementing quantum 
parameter estimation must also develop a constructive 
procedure for doing the estimation. 

In this section, we consider two methods for estimat- 
ing the strength of the magnetic field B based on the 
stochastic measurement record ^(o,t)- In both cases, we 
extend the quantum filters developed in the previous sec- 
tion to account for our uncertainty in 5, which in turn 
results in new filters capable of estimating B. We note 
that our approach is similar in spirit to Refs. |5l|6]. 



A. Quantum Particle Filter 



distribution. This suggests we write the total conditional 
density matrix in the ensemble form 



J dBpt{B)\B){B\®pf^ (54) 



where Pt{B) = P(5|Z(o,t)) is precisely the conditional 
probability density for B. 

While one could attempt to update this state via the 
quantum filter, doing so is entirely impractical, as one 
can not represent an arbitrary distribution for pt{B) with 
finite resources. Instead, one approximates the distribu- 
tion with a weighted set of point masses or particles: 



N 



MB)^Y.Pt^^iB-B,). 



(55) 



The approximation can be made arbitrarily accurate in 
the limit of N ^ oo. Plugging this distribution into the 



ensemble density matrix form of Eq. 54 gives 



N 



(56) 



The technique of quantum particle filtering, as devel- 
oped in Ref. [15 , leverages the fact that the quantum 
filtering equations already provide a means for estimat- 
ing the state of a quantum system conditioned on the 
measurement record. If we place the magnetic field pa- 
rameter on the same footing as the quantum state, we 
can simply apply the quantum filtering results we already 
derived. Indeed, by embedding the magnetic field param- 
eter as a diagonal operator in an auxiliary Hilbert space, 
the quantum filter still gives the best estimate of both 
system and auxiliary space operators. We accomplish 
this by promoting the magnetic field parameter to the 
diagonal operator 



/ 



B^B= B\B){B\dB eHb, 



(52) 



where Hb the new auxiliary Hilbert space with basis 
states satisfying B\B) = B\B) smd {B\B') = S{B - B'). 
All atomic operators and states, which are associated 
with the atomic Hilbert space TYa, act as the identity 
on TYb, e.g. Fz ^ I ^ Fz. The only operator which joins 
the two spaces is the magnetic dipole Hamiltonian, which 
is now given by 



Hi 



-hjB (g) Fy 



(53) 



The derivation of the quantum filtering equation is essen- 
tially unchanged, provided one replaces atomic operators 
with the appropriate forms for the joint space Hb ^ ^a- 
For parameter estimation, the adjoint form is the more 
convenient version of the quantum filter. Since B corre- 
sponds to a classical parameter, we require the marginal 
density matrix {pB)t = ^^HAlPt] be diagonal in the ba- 
sis of B^ so that it corresponds to a classical probability 



Each of the N triples {p[^\ Bi^ Pt^'^} is called a quantum 
particle. Intuitively, the particle filter works by discretiz- 
ing the parameter space and then evolving an ensemble 
of quantum systems according to the exact dynamics for 
each parameter value. The filtering equations below per- 
form Bayesian inference on this ensemble, updating the 
relative probabilities of particular parameter values given 
the measurement record. 

The quantum particle filter for the double-pass system 
with unknown B is found by plugging the discretized 
ensemble pf into the extended double-pass filter. After 
a little manipulation, one finds 



dp['^ = 2Vm{Tt[Fzp[^'^] 



N 



(57a) 



dp 



^MV[Fz]p[^'^dt + KV[Fz]p[^'^dt (57b) 



(yMM[Fz 



Pi 



i^[Fy,p[''^^])dWt 



N 



dWt 



= dZt-2VM^p['^TT[Fzp[^'^]dt (57c) 



where the prior distribution po{B) is used to determine 
the initial parameter weights, Pq'\ and values, Bi. All 

initial quantum states, p^ ' , are taken to be the spin 
coherent state pointing along -\-x. 

Looking at the structure of this filter, we see that each 
particle's quantum state pi evolves under the standard 
quantum filter we would use if we knew B = B^. The 
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innovations process dWt serves to couple the different 
particles, since it depends on the ensemble average of the 
measurement observable. 

An estimate of the magnetic field strength is then con- 
structed from the approximate density in Eq. ([55| , either 
taking the most probable B value, corresponding to the 

(i) 

largest pl or calculating the expected value of B 



N 



(58) 



i=l 



For the latter estimate, the uncertainty is given by 



= ((B^) _ 



1/2 



(59) 



One important feature of particle filters which we will 
mention, but not utilize in this paper, is the potential 
for resampling. As presented, the initial particle ensem- 
ble {pq*\ 5^, Pq^*'^} fixes the possible values of Bi at the 
outset. If there are no initial points near the true value 
of 5, the filter will by construction have difficulty find- 
ing this value. Even if one of the particles is associated 
with the true value of 5, the filter still integrates the 
dynamics for a potentially large number of low weight 
particles which contribute little to the estimate. A way 
around this issue is to resample and create new particles 
during the course of integration, either drawing from the 
discrete distribution {p^^} or using finer and finer reg- 
ular grids to hone in on important regions of parameter 
space. But in both cases, it is not clear how to resample 
the associated quantum states pl *Mn order to improve 
the estimate. Rather than considering such options here, 
we refer the reader to [15 , 31 for more discussion. 

Finally, we note that the computational resources re- 
quired for particle filtering are demanding, generally 
requiring space and time which scale as A^(dim7YA)^, 
though we can get down to NidimHA) if the atomic 
state remains pure. Moreover, the particle filter is biased 
for any finite number of particles, although the variance 
of estimates of random variables using {p^*^} converges 
as N~'^ [31 . Although this will not prohibit us from ex- 
ploring the magnetic field sensitivity scaling of systems 
with collective spins up to F = 140, using a large par- 
ticle set is entirely impractical for real systems with F 
orders of magnitude larger. For such systems, we need a 
low-dimensional state representation whose size does not 
scale with F, after which the particle filter becomes a 
more promising technique for parameter estimation. We 
explore low dimensional representations of the double- 
pass system in the following subsection. 



B. Quantum Kalman Filter 

Rather than constructing a magnetic field estimator 
from the exact quantum dynamics, one could instead first 
focus on deriving an approximate filter for the atomic 
state, which is then a starting point for the magnetic field 
estimator. Indeed, previous work in precision magnetom- 
etry via continuous measurement [5^ has taken this route 
by constructing a quantum Kalman filter to describe the 
atomic dynamics. Such a filter leverages the fact that for 
an initially spin polarized state of many atoms (say along 
a first order Holstein-Primakoff expansion [19 lin- 
earizes the small-angle dynamics in terms of a Gaussian 
state characterized by the means tt^ [Fz], 7rt[Fy] and the 
covariances AF^^, AFy , AF^Fy. Just as in classical filter- 
ing theory, the conditional state for a linear system with 
Gaussian noise is itself described by a Gaussian distri- 
bution and therefore only requires filtering equations for 
the means and a deterministic equation for the variances 
[T7[ HH]. For the case of magnetometry, the number of 
these parameters is independent of the number of atoms 
in the atomic ensemble. We will also find that within this 
approximation, we can again embed B as an unknown 
state parameter and find a corresponding Kalman filter 
appropriate for estimating its value. 

However, applying the small-angle and Gaussian ap- 
proximations in the quantum case is usually done in an 
ad-hoc fashion, especially in light of the recent introduc- 
tion of projection filtering into the quantum filtering set- 
ting [20l[2T]. In this framework, one selects a convenient 
manifold of states whose parameterization refiects the ap- 
proximations to enforce. At each point in this manifold, 
the exact differential dynamics induced on these states is 
orthogonally projected back into the chosen family. For 
our purposes, this means projecting the filter in Eq. (46) 



onto a manifold of Gaussian spin states. Although the 
resulting equations are not substantively different than 
those derived less carefully, we believe the potential ap- 
plication of projection filtering in deriving other approxi- 
mate filters and master equations warrants the following 
exposition. 



1. Projection Filter Overview 

Abstractly, projection filtering proceeds as follows. We 
assume we already have a dynamical equation, such as 
Eq. (46), for a given manifold of states, such as pure 
states. For convenience, let these dynamics be repre- 
sented as 



(60) 



where M is the generator of dynamics. Now select the 
desired family of "approximating" states which are a 
submanifold of the exact states. We assume this fam- 
ily is parameterized by a finite number of parameters 
xi,X2,...,Xn and we denote states in this family as 
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X2, . . . , At every point in this manifold, the tan- 
gent space is spanned by the tangent vectors 



dxi 



(61) 



Loosely speaking, these tangent vectors tell us how dif- 
ferential changes in the parameters move us through the 
corresponding submanifold of ^2, . . . , x^) states in 
the space of pure states. This is particularly useful, as 
the action of the generator A/'[|xi,X2, . . . ,^n)] does not 
necessarily result in a state within the family. But by 
projecting the dynamics onto the tangent space, we can 
find a filter, called the projection filter, which constrains 
evolution within the chosen submanifold. Explicitly, this 
projection is written as 



fa,A/'[|xi,X2,...,Xn)]) ^ 



(62) 



where in this pure state formulation, the inner product 
is the standard Hilbert space inner product. 



2. Gaussian State Family and Tangent Vectors 

For our double-pass magnetometer, we begin by intro- 
ducing the two-parameter family of Gaussian states 

= f,,%|F,+F,), (63) 

where \F^^Fx) is the spin coherent state pointing along 
+x, S^^ is a spin squeezing operator with squeezing 
parameter and Yq^ is a rotation about the ?/-axis by an- 
gle Ot. Intuitively, the squeezing along z generated by S^^ 
corresponds to the squeezing induced by measuring F^. 
The rotation via Yq^ then accounts for both the random 
evolution due to the measurement as well as any rotation 
induced by the magnetic field. The tangent vectors for 
these states are 



d\eu^t) 

dOt 

-iF^YeAAF.^F,) 
d\Ou^t) 



(64) 



= f,,4(-2z(F,Fy + FyF,))|F,+F,). (65) 

In calculating the normalization of these tangent vec- 
tors, we encounter terms such as 

{ve,,ve,) = {F, +F,\sIf^S^, \F, +F,). (66) 

More generally, almost all inner-products needed for the 
projection filter will be of the form 



Here, g and / are polynomial functions of their argu- 
ments. Since Yq^ is a rotation, we can exactly evaluate 



Ylf{F^,F^,F,)Ye, = f{YlF^Yo„YlFyYe,,YlF,YoJ, 
where 



(67) 



Yo\F^Yg^ = FM = cos Ot + F, sin Ot 
^Ot^y^Ot = Fy 

YlP^Ye^ = F,(^,) = ^zCos^,-FxSin^^ 
This leaves us with expectations of the form 



(F, +F,\Slg{F^, Fy, F^)fiFM, Fy, ^z(^*))% 1^^, +Fcc) 

(68) 

where gx f will just be linear combinations of powers and 
products of Fx,Fy,Fz. Unfortunately, we cannot evalu- 
ate this expectation for arbitrary ^t- However, for small 
the state which we are taking expectations with re- 
spect to is the "squeezed vacuum" in our preferred basis, 
e.g. it is the state |F, +Fc) pointing in the same direc- 
tion, but with squeezed uncertainty in Fz and increased 
uncertainty in Fy. 

For large F, angular momentum expectations of such 
a state are extremely well described by the Holstein- 
Primakoff approximation to lowest order [19] 



2Fa 
2Fa^ 



(69) 



where F±^x = Fy ± iF^^ and a^a^ are bosonic creation 
and annihilation operators. We then write our state as 
|F, +Faj) = |0), which is the vacuum in the Holstein- 
Primakoff representation. Under this approximation, we 
can use the relations 



= F 



e4^«*(a + at) 



/2F 



* (a-a'l') 



(70) 
(71) 
(72) 



to evaluate the expectation in Eq. (68). In light of this 



approximation, the tangent vector overlaps are readily 
shown to be 



: 8F^ 

0, 



(73) 

(74) 
(75) 



where the last result indicates the tangent vectors are 
orthogonal as desired. 
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3. Orthogonal Projection of Double-pass Filter 

Before performing orthogonal projection of the dynam- 
ics onto the tangent space, we must first convert the 
filtering equation from Ito to Stratonovich form. As is 
discussed in Ref. [27 , the Ito chain rule is incompatible 
with the differential geometry picture of projecting onto 
the tangent space. Fortunately, Stratonovich stochas- 
tic integrals follow the standard chain rule and are thus 
amenable to projection filtering methods. Following the 
derivation in Appendix [C) we find that the Stratonovich 
SSE is given by 



Chugging through the remaining terms in a similar fash- 
ion, we arrive at the full projection filter equations 



dOt -fBdt- 



and 



sin Otdt + 2FVkM sin Otdt 



Me-^^^* cos6>t + 



odWt 



(82) 



d^t = ^e-'^^^ cos^Mt. 



(83) 



(F, - {F^),f - {AF^), 



^2iVKM{F,)tFy + iVKM{F,Fy)t \^)tdt 



(Vm(F, - (F,)t) + zV^Fy) \^)t o dWu 



(76) 



where {AF^)t = (F^) - {F,)\ 

In order to find the projection filter, we compare the 
general projection formula in Eq. ( [62| ) to the general dy- 
namical equation for states in our chosen family, given 

by 



(77) 



Using the orthogonality of the tangent vectors, the gen- 
eral forms for d^t and d9t are 



d0t 
d^t 



F 



8F2 



{vi„d\iJtmt,Ot]), 



(78) 
(79) 
under the 



where d\tljt)[£,t^^t] is the evolution of \£,t^^t) 
Stratonovich filter of Eq. (76). 

As an example calculation using these methods, con- 
sider projecting the dynamics generated by the magnetic 
field term. Its contribution cq to the Ot dynamics is given 
by 



ce 



F 
F 



{v0,,-ijBFy\et,^t)dt) 



(0|5^^f;/y2f,,%|0)dt (80) 



= jB{0\{a^a'^f\0)dt 
jBdt. 

Similarly, the contribution to is 

= ^{^^t^-nBFy\Ot,it)dt) 
-fB 



^{0\S^:{F^Fy^FyF^)YlFyYe,S^jO)dt 



= 0. 



-a^a"^ - a"^ a - a"^ \0)dt 



(81) 



Converting back to Ito form using Eq. ( |C3[ ), we have 
dOt 



M 



Bj - — e-^^^^* sin(2(9t) + 2FVKMsinet 



Me-^^^' cos Of 



dWt 



dit 



M 



e"^^^* cos^ Otdt, 



dt 
(84a) 
(84b) 



where the innovations are now in terms of the approxi- 
mation of (Fz)t within the Gaussian family: 



(85) 



dWt = dZt-2VM{F^)tdt 
= dZt^2FVMsin0tdt. 

4' Small- angle Kalman Filter 



We see that the projected filter in Eq. (84) is actually 
more general than the filters usually derived for the mag- 
netometry problem, which do not distinguish the Gaus- 
sian and small- angle approximations. That is, the family 
of states in Eq. ( 63 ) and the approximations considered in 
the above derivation only enforce the Gaussian state as- 
sumption through the Holstein-Primakoff approximation. 
We can separately apply the small-angle approximation 
to recover an equation appropriate for the Kalman filter. 
In this limit, the equation for completely decouples 
and has a closed form solution 



6 = ^ln[l + 2™]. 



(86) 



Taking the small- angle approximation for 6t and plugging 
in the explicit form of gives 



dOt 



Bj + (2FVkM 



M 



2(1 + 2FMtY 



1 + 2FMt 



dt 



dWu (87) 



which is linear in the remaining state parameter 9t. 

While we could consider the Kalman filter for the quan- 
tum state alone, we can just as easily account for our 
uncertainty in B at the same time. That is, if we now 
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embed 5 as a state variable, setting Xf = [Of B]^ ^ the 
dynamics can be written in a linear form as 



dXt 
dZt 

A 



AXtdt + BdWt 
CXtdt + DdWt 



f2Fy/KM 



M 







2(1+2FM^)2 



B = 

C = 
D = 



7m 

l-\-2FMt 





O) 

1. 



(89) 
(90) 

(91) 

(92) 
(93) 



Equations (^88l and (89) are precisely a classical linear 



system/observation pair, in which the same white noise 
process (the innovations) drives both the system and ob- 
servation processes. The estimate Xt = E[Xt|Z(o,t)] ad- 
mits a Kalman filter solution [33 , given by 

dXt = AXtdt^{B^VC'^)dWt (94) 
V = AV ^VA"^ ^ BB"^ - {B ^VC'^){B ^VC'^y 

where V is the covariance matrix 

V = E[(X - E[X])(X - E[X])^] 

MBkf ABlf J 

and 

2FVM0tdt 



(95) 
(96) 



dWt = dZt 



(97) 



is the innovations constructed from the current Ot esti- 
mate in the small- angle approximation. 

Looking at the explicit system of equations for the vari- 
ances, we have 



dt 



-MAOl 



+ 4F + ^F'^Mt 
(1 + 2FMtY 



d{ABl,) 



kfJ _ 



dt 

djAOtBkf) 
dt 



-^2-fAOtBkf 
-4F^M{A0tBkff 



AF'^AOl 
(98) 
(99) 



M 



kf 



2(l + 2FMt)2 
(1 + 4F + SF^Mt^ 

8F2(1 + 2FMtfA0l) AOtBkf 



(100) 



which are completely independent of the second-pass cou- 
pling strength K. That is, within the small-angle and 
Gaussian approximations, the double-pass system has no 
improvement in sensitivity and gives rise to the same 
F~^ uncertainty scaling found previously for single-pass 
systems [5]. Perhaps this is unsurprising, as we at- 
tempted to find a linear description of an essentially 
non-linear affect. Indeed, the numeric simulations in the 



next section suggest the single-mode Gaussian approx- 
imation breaks down just as the double-pass filter be- 
gins to show improved sensitivity to the magnetic field 
parameter. Finally, we note that we have also derived 
a filtering equation which retains the next term in the 
Holstein-Primakoff expansion, but whose K dependence 
nonetheless shows a negligible change relative to the low- 
est order expansion. 



V. SIMULATIONS 

Given the absence of an analytic improvement in the 
sensitivity of the quantum Kalman filter, we turn to nu- 
merical simulations of the quantum particle filter in or- 
der to gauge the potential of the double-pass system for 
magnetometry. First recall how the filter would be used 
in an actual experiment. Continuous measurements of 
the atomic cloud Larmor precessing under a particular, 
albeit unknown, magnetic field B would give rise to the 
observations process ^(o,^)- This would then be fed into a 
classical computer to propagate the quantum particle fil- 
tering equations given in (57). The computer would then 
use the quantum particle set to provide the estimate Bpf 



and uncertainty AB'^^. 

In order to simulate such an experiment, we can gen- 
erate the stochastic measurement record ^(o,t) using the 
quantum filter for the double-pass system given in Eq. 
( [46| ), evolved with a known magnetic field B. Since the 
system is driven by the white noise process dWt^ the fil- 
tering equations may be integrated by the same integra- 
tor previously used to approximate the quantum Cramer- 
Rao bound. The measurements generated by these tra- 
jectories are equivalent to what the quantum particle fil- 
ter would receive in an experiment, which means they 
can then be fed into the same particle filtering code to 
simulate an estimate oi B. In order to compare perfor- 
mance, we actually simulate two systems in parallel, one 
representing the double-pass system and the other, with 
= 0, representing a single-pass system. Both utilize 
the same noise realizations on an individual trajectory. 

As is common when considering detection limits, we 
focus on the case of 5 = 0. Although an unbiased esti- 
mator would assume no prior knowledge of the magnetic 
field value, such an approach is impractical for the par- 
ticle filter, which would fail in approximating such large 
uncertainty with a finite number of particles. As such, we 
take the initial distribution of B values for the quantum 
particle set to be Gaussian 



Pq{B]I1b,CFb) 



1 



: exp(- 



2al 



(101) 



with mean /j^b = ^ and variance a% = lOi^^, where we 
again set 7 = 1 and again define all parameters in units 
of u. For a set of N particles, the particle magnetic field 
values {Bi} are drawn from the initial distribution, with 

weights = 1/N. The initial quantum state for all 
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FIG. 3: (color online) (a) Estimator uncertainties as a function of F averaged over 100 trajectories with M — lOz^, K — 0.0006z^, 
B — ^ and ts — 0.1z^~^. The initial N — 1000 particle set was drawn from a Gaussian distribution with mean zero and variance 
lOz^^, which was also the same initial uncertainty in the Kalman filter ABkf. A power- law fit to the particle filter (PF) scalings 
shows a single-pass scaling of and a double-pass scaling of Also shown are the quantum Cramer-Rao (QCR) 

bounds previously simulated for Figure [5] The inset shows the sample estimator deviation Spf for the same simulations, (b) 
Average effective particle number fraction Nefi{ts)/N as a function of F for the data in (a). 



particles is set to the spin-coherent state along +x,i.e. 

Figure [3]^ a) shows, with solid lines, the average parti- 
cle filter uncertainty ABpf as a function of F, averaged 
over 100 measurement realizations using N = 1000 par- 
ticles in each run of the filter. The error bars represent 
the the deviation in the simulated uncertainties over the 
100 runs. As was the case for the Fisher information 
calculations, we observe an improved sensitivity scaling 
for the double-pass system, albeit with increased fluctua- 
tions in the individual run uncertainty ^Bpf. Power-law 
least-squares fits of the average give a single-pass uncer- 
tainty scaling F~^-^^ and a double-pass scaling of F~^-^^ 
which are consistent with the quantum Cramer-Rao scal- 
ings in figure |2j Also shown is the analytic single-pass 
uncertainty scaling given by numerical integration of the 
Kalman covariance matrix via Eq. (98). We see that this 



agrees very well with the single-pass particle filter scal- 
ing and since it is consistent with previous Kalman filters 
used for magnetometry [5 , suggests the double-pass scal- 
ing does indicate improved sensitivity. 

Of course, these statements are not without caveats. 
The dashed lines in the plot correspond to the numer- 
ically computed quantum Cramer-Rao bound, which is 
clearly below the estimates of all the filters. This might 
mean that the continuous-measurement which gives rise 
to the numerical bound via Eq. (40) is simply not satu- 



rated by the corresponding estimator for that continuous- 
measurement. Unfortunately, the above data took a 
week to generate on a quad-core workstation, indicating 
the technical challenges already present in simulating an 
= 1000 quantum particle set for the depicted range 



of F limits the quality of the statistics. As previously 
mentioned, the particle filter approximation is inherently 
biased, with the variance of estimates converging as N~^. 
The inset in figure [sja) shows the sample estimator de- 
viation Spf^ which is the deviation in the actual perfor- 
mance error of the particle filter on each individual run, 
i.e. Bpf — B where the true B = 0. In other words, ^Bpf 
is the uncertainty calculated for an individual trajectory 
from the particle distribution {Pt*'^}, which is averaged 

over many trajectories to get ^Bpf. However, an in- 
dividual run of the particle filter also gives an estimate 
Bpf of the true magnetic field B. Since we know that 
the measurements were generated from a system evolved 
with 5 = 0, we can calculate the deviation in the ac- 
tual estimates Bpf. If the particle filter were unbiased, 
we would expect this sample deviation to equal the av- 
erage particle filter deviation, i.e. Spf = ^Bpf. Instead, 
the sample deviation dwarfs the average estimator uncer- 
tainty, indicating that the particle filter bias dominates. 

We believe this is primarily a technical issue due to 
the dwindling contribution of initial particles to the final 
estimate. Given that the potential particle magnetic field 
values {Bi} are fixed at the outset, particles far from the 
true value will contribute only marginally to the final 
estimate. That is, for Bi far from the true 5 = 0, pi^ 
will be relatively small. This behavior is well known in 
classical particle filtering [34^ and is characterized by the 
effective sample size 



iVeff(t) 



(102) 
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Figure [3|^b) shows the average effective particle number 
fraction A^eff(^s)/^ and its deviation as a function of N 
for the same trajectories used in part (a) of the figure. 
By the end of integration, we see that both single and 
double-pass systems have less than 100 of the initial 1000 
particles significantly contributing to the estimate. The 
fact that the double-pass system shows a smaller effective 
particle number is consistent with the increased sample 
estimator deviation Spf. In order to decrease the bias, 
one would need to include more particles or develop a re- 
sampling technique. Nonetheless, we still believe the sim- 
ulations depicted in figure [sja) suggest that the double- 
pass system shows in improved uncertainty scaling. 

Numerical simulation also provides insight into how 
the Gaussian state assumption of the Kalman filter ap- 
plies in the double-pass case. Figure [4] shows quantum 
states evolved under two different noise realizations with 
B = 0,M = lOu, K = 0.0006iy. Both states were ini- 
tially spin-polarized along -\-x and evolved until time 
ts = 0.1u~^ under the fuh double-pass SSE in Eq. (46). 
The Q-function shown is defined as 



Q{e,^,t) = \{0,^\i^t)\' 



(103) 



where the spin-coherent state is the -\-F eigenstate 

of the spin-operator Fx sin cos (j^^Fy sin sin (j)^F^ cos 0. 
Although one example shows a Gaussian squeezed spin 
state, the other shows a state with a bimodal Gaussian 
distribution, which we suspect is a consequence of the 
complicated nonlinear dynamics of Eq. (|46| — the pres- 
ence of operators that are nonlinear in the spin observ- 
ables, such as F^Fy^ do not preserve Gaussian states. 
Therefore, a filter which confines the dynamics to the 
Gaussian family in Eq. (63) is only likely to remain valid 



for short time, and helps explain why the Kalman filter 
fails to find a difference between the single and double- 
pass setup. These plots suggests a family of bimodal 
Gaussian states might result in a useful projection fil- 
ter; we have been unable to find a parameterization of 
such a family which admits an analytic derivation of a 
projection filter. 



VI. CONCLUSIONS 

We have explored the use of double-pass continuous 
measurement for precision magnetometry. Our primary 
result involves numerical simulations of the quantum 
Cramer-Rao bound which indicate that a double-pass 
system shows an improved magnetic field uncertainty 
scaling with atom number over a comparable single-pass 
system, albeit only for particular choices of coupling 
strengths relative to the collective spin size. This is in 
contrast to quantum information theoretic bounds which 
suggest that the Heisenberg limit bounds the uncertainty 
scaling for both a single and double-pass system. Clearly, 
future work aimed at reconciling these results is neces- 
sary, particularly deriving analytic quantum Cramer- Rao 
bounds for unbounded ancilla systems. However, at a 
heuristic level, we see that the double-pass provides a 
positive feedback effect which amplifies the amount of 
Larmor precession due to the magnetic field. This of 
course also amplifies the initial spin-projection uncer- 
tainty, but the continuous measurement allows us to over- 
come this uncertainty by actually learning the state of 
the collective spin. It is this interplay between amplifica- 
tion and measurement-induced squeezing which requires 
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a fine-tuning of parameter strengtiis in order to find an 
improved Cramer-Rao bound. 

We tiave also explored estimators intended to achieve 
the uncertainty scaling seen in numerical simulations. 
Taking a brute force approach, quantum particle filters 
show evidence of the improved double-pass scaling, al- 
though the results suffer from limited statistics which can 
not be significantly improved with current computational 
power and methods. We are likewise limited to study- 
ing relatively small atomic ensemble sizes (small collec- 
tive spin sizes), though we expect the improved scaling 
to extend to larger ensembles. More practical quantum 
Kalman filters show no improved sensitivity, which is con- 
sistent with an observed breakdown in the Gaussian state 
assumption used to derive them. However, the general 
projecting filtering technique used in the Kalman filter 
derivation provides an avenue for deriving more appro- 
priate filters which might prove more tractable for practi- 
cal magnetic field estimation. More generally, we believe 
effective nonlinear interactions may prove an important 
tool in precision measurement. 
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useful disagreements and discussions. This work was 
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APPENDIX A: QUANTUM ITO PRODUCT 
RULE 

When manipulating products of variables governed by 
stochastic differential equations, one must be careful to 
use the Ito product rule [35 , which generalizes for the 
quantum case j29l|36j as follows. Consider two operators 
governed by the quantum stochastic differential equa- 
tions 

dXt = AdSt + BdSl + Cdt (Al) 
dYt = bdSt + EdSl + Fdt (A2) 

then the product X^F^ is governed by 

d{XtYt) = Xt[dYt) + {dXt)Yt + {dXt){dYt) 

= {XtD + AYt)dSt + {XtE + BYt)dSl 

+ {XtF + CYt + AE)dt, (A3) 

where we have used the fundamental Ito products — 
dStdSl = dt and dSldSt = dS^ = {dS})'^ = 0. Some 
find it useful to heuristically identify the white noise dif- 
ferentials dSt^dSl as order ^/dtt^ so that any consistent 
chain rule requires keeping terms to second order in dif- 
ferentials or equivalent ly, first order in dt. 



APPENDIX B: STOCHASTIC SCHRODINGER 
EQUATION 

Lacking any extra sources of decoherence, pure states 
remain pure under the dynamics described by the quan- 
tum filtering equation. As such, it is often convenient for 
analysis and simulation to have a pure state description of 
the dynamics in terms of a stochastic Schro dinger equa- 
tion (SSE). Although these have previously appeared in 
the literature, for completeness we briefly derive the SSE 
for the general adjoint filter 

dpt = -i[H,pt]dt^ (iptL^ - ^L^Lpt - \ptL^i^ dt 

+ (Lpt + ptD - Tr [(L + D)pt]pt) dWt. (Bl) 

We begin by writing 

d\^)t = A\^)tdt^B\^)tdWt (B2) 
d{^\t = {^\tAUt + {^\tBUWt (B3) 

From the Ito rules, we have 

d{pt) = d{\^){^\t) 

= |^/')td((^/;|t)+d(|^/;)t)(V'|t + d(|^/'t))d((^/'k) (B4) 
= [Apt + ptA^)dt + {Bpt + ptB^)dWt + BptB^dt 

Comparing the coefficients to the quantum filtering equa- 
tion, we read off 

B = L-{L) (B5) 
fit = it _ ^j;^t) (B6) 

so that 

BptB^ = LptL^ - {L^)Lpt - {L)ptL^ + {L){L^)pt (B7) 
We try setting 

A = -i(LtL-2{Lt)L+(L)(Lt)) (B8) 

which means that 

Apt + ptA^ = -\ {L^L - 2{L^)L + {L){L^)) pt 

-pt\{L^L-2{L)L^ + {L){L^)) (B9) 

= -h^Lpt-^PtL^L + {L^)Lpt 

+ PtL^L)-{L){L^)pt (BIO) 

and therefore 

Apt + ptA^ + BptB^ = -^L^Lpt - \ptL^L 

+ {L^)Lpt+ptL^L) 

-{L){L^)pt + LptL^ 

-{L^)Lpt-{L)ptL^ + {L){L^)pt 

(Bll) 

= LptL^ - h^Lpt - \ptL^L 

(B12) 
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which is the deterministic part of the quantum filtering 
equation as desired. 



APPENDIX C: CONVERTING BETWEEN ITO 
AND STRATONOVICH SDE 

Since white noise is non- different iable, there are a vari- 
ety of ways to give meaning to the differential dWt when 
integrated against some other stochastic process. The 
most statistically pleasing is the Ito integral, represented 
by the SDE 



dxt = Xt)dt + Xt)dWt 



(Cl) 



and which satisfies (iE[x^] = a(t^Xt)dt^ but requires us- 
ing the Ito chain rule. Conversely, a Stratonivich SDE, 
written 



dxt = a(t^ Xt)dt + Xt) o dWt^ 



(C2) 



transforms using the normal chain rule but has a non- 
trivial expectation. Fortunately, it is straightforward to 
convert between the two forms, which share the same 
stochastic coefficient and whose deterministic coefficients 
are related via 



1 

a^(t,xt) =a\t,xt) - - ^h^(t,xt) 



dW{t,Xt) 



k=l 



dx^ 



(C3) 



where the superscripts denote the j-th or k-th entry in 
the corresponding vector. We recommend the interested 
reader consult [35] for a more complete discussion of 
the various ways of constructing stochastic integrals and 
defining white noise increments. 



For the double-pass Ito SSE in Eq. (46), we begin 



the conversion by noting that states with entirely real 
amplitudes form an invariant set and therefore write 
\i^)t = Xlm=-F^n^)- stochastic coefficient is then 



b{t,Xt) = VM mxr\m)-y/M ^ n{xt''fxr\ 

m= — F m,n— — F 

1 ^ 

+ -V^ [V{F - m){F + m + l)xr\m + 1) 



ra^-F 



+ m){F - m + l)xr\m - 1) 



(C4) 



which has as its j-th entry 



V{t,Xt) = VM{j- J2 nixtn^)xt^ 



n=-F 



^/{F-j + l){F + j)xt 



i-1 



(C5) 



-^{F + j + l){F-j)xi+^ 
The derivative with respect to Xt^ is then 

db^{t,xt) 



dxk 



Y n{xt''f)Sjk-^2kxt^Xt' 

K 



n^-F 



V(i^-i + l)(F+j)%-i),fc 



SO that the sum in Eq. (C3) is 



k=-F 



dXf^ 



= s/M{j- n{xt^f)V{t,Xt)-2s/MY,kxtH\t,Xt) 



Xt' 



(C6) 



y^F-3 + l){F + j)V-\t,Xt) - V{F + j + l){F-j)V+\t,Xt) 



This suggests an equivalent operator form 



MiF,-{F,)t)+i^FyY -2VM{F,{^/MiF,-{F,))+i^/KFy))t = \Vm{F, - {F,)t) + i^Fyf (C7) 



-2M{AF^)t - 2i^/KM{F^Fy) 



y/i> 



where {AF^)t = (F^) - {F^f, so the Stratonovich SSE is 



-i'jBFy - M 



(F, - {F,)tf - {AF^)t \ - ^^^4 + 2iVKM{F,)tFy + i^KM{F^Fy)t 



\iP)tdt 



(yM{F, - {F,))t + iVi^Fy) \i^)t o dWt. 



(C8) 



18 



[1] D. Budker, W. Gawlik, D. F. Kimball, S. M. Rochester, 
V. V. Yashchuk, and A. Weis, Rev. Mod. Phys. 74, 1153 
(2002). 

[2] C. Helstrom, Quantum Detection and Estimation Theory 

(Academic Press, New York, 1976). 
[3] S. L. Braunstein and C. M. Caves, Phys. Rev. Lett. 72, 

3439 (1994). 

[4] I. K. Kominis, T. W. Kornack, J. C. Allred, and M. V. 

Romalis, Nature 422, 596 (2003). 
[5] J. Geremia, J. K. Stockton, A. C. Doherty, and 

H. Mabuchi, Phys. Rev. Lett. 91, 250801 (2003). 
[6] K. M0lmer and L. B. Madsen, Phys. Rev. A 70, 052102 

(2004). 

[7] A. Holevo, Probabilistic and Statistical Aspects of Quan- 
tum Theory (North-Holland, Amsterdam, 1982). 
[8] S. L. Braunstein, C. M. Caves, and G. Milburn, Ann. 

Phys. (N.Y.) 247, 135 (1996). 
[9] S. Boixo, S. T. Flammia, CM. Caves, and J. Geremia, 
Phys. Rev. Lett. 98, 090401 (2007). 

[10] D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. 
Heinzen, Phys. Rev. A 50, 67 (1994). 

[11] A. M. Rey, L. Jiang, and M. D. Lukin, Phys. Rev. A 76, 
053617 (pages 8) (2007). 

[12] J. F. Sherson and K. M0lmer, Phys. Rev. Lett. 97, 
143602 (pages 4) (2006). 

[13] G. Sarma, A. Silberfarb, and H. Mabuchi, Phys. Rev. A 
78, 025801 (2008). 

[14] C. A. Muschik, K. Hammerer, E. S. Polzik, and J. L 
Cirac, Phys. Rev. A 73, 062329 (2006). 

[15] B. A. Chase and J. M. Geremia, Phys. Rev. A 79, 022314 
(2009). 

[16] V. Belavkin, Rep. Math. Phys. 43, 405 (1999). 
[17] R. Kalman, Trans. ASME Ser. D. J. Basic Eng. 82, 35 
(1960). 

[18] R. E. Kalman and R. Bucy, Trans. ASME Ser. D. J. Basic 

Eng. 83, 95 (1961). 
[19] F. Holstein and H. Primakoff, Phys. Rev. 58, 1048 

(1940). 

[20] R. van Handel and H. Mabuchi, J. Opt. B: Quantum 

Semiclass. Opt. 7, S226 (2005). 
[21] H. Mabuchi, Phys. Rev. A 78, 015801 (2008). 
[22] G. A. Smith, S. Chaudhury, A. Silberfarb, L H. Deutsch, 

and P. S. Jessen, Phys. Rev. Lett. 93, 163602 (2004). 
[23] J. M. Geremia, J. K. Stockton, and H. Mabuchi, Phys. 

Rev. A 73, 042112 (pages 14) (2006). 
[24] L. Accardi, Y. Lu, and L Volovich, Quantum Theory and 

Its Stochastic Limit (Springer, 2002). 
[25] L. Accardi, A. Frigerio, and Y. Lu, Comm. Math. Phys. 

131, 537 (1990). 
[26] J. Cough, Commun. Math. Phys. 254, 489 (2005). 
[27] R. van Handel, J. K. Stockton, and H. Mabuchi, J. Opt. 

B: Quantum Semiclass. Opt. 7, S179 (2005). 
[28] L. Bouten, J. Stockton, G. Sarma, and H. Mabuchi, Phys. 

Rev. A 75, 052111 (2007). 
[29] L. Bouten, R. van Handel, and M. James, SIAM J. Con- 
trol Optim. 46, 2199 (2007). 
[30] P. E. Kloeden and E. Platen, Numerical Solution of 

Stochastic Differential Equations (Springer, 1999), see 

pages 503-506. 

[31] A. Doucet, N. de Freitas, and N. Gordon, eds., Sequential 
Monte Carlo Methods in Practice (Springer, 2001). 



[32] M. Kitagawa and M. Ueda, Phys. Rev. A 47, 5138 (1993). 

[33] R. Lipster and A. Shiryayev, Statistics of Random Pro- 
cessses I: General Theory (Springer- Verlag, 1977). 

[34] M. Arulampalam, S. Maskeh, N. Gordon, and T. Clapp, 
IEEE Transactions on Signal Processing 50, 174 (2002). 

[35] B. 0ksendal, Stochastic Differential Equations: An In- 
troduction with Applications (Springer, 2002). 

[36] K. Parthasarathyy, An Introduction to Quantum Stochas- 
tic Calculus (Birkhauser, 2002). 



